Fit example models

Author

Max Lindmark

Published

2024-04-08

# Load libraries
library(tidyverse)
library(tidylog)
library(sdmTMB)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(modelr)
library(ggstats)
library(ggsidekick); theme_set(theme_sleek())

home <- here::here()
source(paste0(home, "/R/functions/map-plot.R"))

Read data and prediction grid, scales variables

# Read data
d <- readr::read_csv(paste0(home, "/data/clean/larval_size.csv")) %>% 
  drop_na(temp) %>% 
  mutate(temp_sc = (temp - mean(temp, na.rm = TRUE)) / sd(temp, na.rm = TRUE),
         yday_ct = yday - mean(yday),
         year_f = as.factor(year),
         species_f = as.factor(species),
         year_ct = year - median(year))

# Trim some outliers
d %>% 
  group_by(species) %>% 
  mutate(sd = sd(length_mm),
         mean = mean(length_mm)) %>% 
  ungroup() %>% 
  mutate(outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
                          "No", "Yes")) %>% 
  ggplot(aes(year, length_mm, color = outlier)) + 
  facet_wrap(~species, scales = "free") + 
  geom_point()

d %>% 
  group_by(species) %>% 
  mutate(sd = sd(length_mm),
         mean = mean(length_mm)) %>% 
  ungroup() %>% 
  mutate(outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
                          "No", "Yes")) %>% 
  ggplot(aes(length_mm, fill = outlier)) + 
  facet_wrap(~species, scales = "free") + 
  geom_bar()

sort(unique(d$species))
 [1] "Agonus cataphractus"      "Ammodytidae"             
 [3] "Anguilla anguilla"        "Aphia minuta"            
 [5] "Argentina silus"          "Argentina spp"           
 [7] "Chirolophis  ascanii"     "Clupea harengus"         
 [9] "Crystallogobius linearis" "Enchelyopus cimbrius"    
[11] "Limanda limanda"          "Maurolicus muelleri"     
[13] "Microstomus kitt"         "Myoxocephalus scorpius"  
[15] "Pholis gunnellus"         "Pomatoschistus sp"       
[17] "Sardina pilchardus"       "Sprattus sprattus"       
[19] "Syngnathus rostellatus"   "Taurulus bubalis"        
d <- d %>%
  group_by(species) %>% 
  mutate(sd = sd(length_mm),
         mean = mean(length_mm)) %>% 
  ungroup() %>% 
  mutate(
    outlier = ifelse(length_mm > mean - 4*sd & length_mm < mean + 4*sd,
                     "No", "Yes"),
    outlier = ifelse(species == "Limanda limanda" & length_mm > 50, "Yes", outlier)) %>%
  filter(outlier == "No")

ggplot(d, aes(length_mm, fill = outlier)) + 
  facet_wrap(~species, scales = "free") + 
  geom_bar()

d %>% 
  group_by(year, species) %>% 
  summarise(n = n()) %>% 
  ggplot(aes(year, n)) + 
  facet_wrap(~species, scales = "free_y", ncol = 3) + 
  geom_bar(stat = "identity")

coul <- brewer.pal(11, "Spectral") 
coul <- colorRampPalette(coul)(length(unique(d$species)))

p1 <- ggplot(d, aes(log(length_mm), fill = species)) + 
  geom_histogram() + 
  scale_fill_manual(values = coul, name = "Species") + 
  coord_cartesian(expand = 0) + 
  labs(y = "Count", x = "Length (mm)", tag = "b)") + 
  theme(legend.text = element_text(face = "italic", size = 7),
        legend.key.size = unit(0.25, "cm"),
        legend.position.inside = c(0.2, 0.69),
        plot.tag = element_text()) +
  guides(fill = guide_legend(position = "inside"))

# Load prediction grid
pred_grid <- readr::read_csv(paste0(home, "/data/clean/pred_grid.csv")) %>% 
  drop_na(temp) %>% 
  filter(year %in% unique(d$year)) %>% 
  mutate(temp_sc = (temp - mean(d$temp, na.rm = TRUE)) / sd(d$temp, na.rm = TRUE),
         year_f = as.factor(year),
         year_ct = 0,
         yday_ct = 0) %>% 
  mutate(keep = ifelse(lon < 10 & lat < 57.15, "N", "Y")) %>% 
  filter(keep == "Y") %>% dplyr::select(-keep)

# Plot temperature in space
p2 <- plot_map + 
  geom_raster(data = pred_grid %>% 
                group_by(X, Y) %>% 
                summarise(mean_temp = mean(temp)),
              aes(X*1000, Y*1000, fill = mean_temp)) + 
  geom_point(data = d %>% distinct(X, Y), aes(X*1000, Y*1000),
             size = 0.5, alpha = 0.3, shape = 21, fill = NA, color = "white") +
  scale_fill_viridis(option = "magma", name = "Temperature (°C)") + 
  theme(legend.position.inside = c(0.25, 0.1),
        legend.direction = "horizontal",
        legend.key.width = unit(0.5, "cm"),
        legend.key.height = unit(0.2, "cm"),
        legend.text = element_text(size = 6),
        legend.title = element_text(size = 8),
        plot.tag = element_text()) +
  labs(tag = "a)") +
  guides(fill = guide_colorbar(position = "inside",
                               title.position = "top",
                               title.hjust = 0.5)) + 
  geom_sf() +
  annotate("text", label = "Sweden", x = xmin2 + 0.95*xrange, y = ymin2 + 0.75*yrange,
           color = "gray50", size = 3) + 
  annotate("text", label = "Norway", x = xmin2 + 0.08*xrange, y = ymin2 + 0.95*yrange,
           color = "gray50", size = 3) +
  annotate("text", label = "Denmark", x = xmin2 + 0.42*xrange, y = ymin2 + 0.45*yrange,
           color = "gray50", size = 3)
  
p2 + p1

ggsave(paste0(home, "/figures/data_map.pdf"), width = 22, height = 12, units = "cm")

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/00-fit-models_cache/html")

Fit alternative models

mesh <- make_mesh(d,
                  xy_cols = c("X", "Y"),
                  cutoff = 5)
    
ggplot() +
  inlabru::gg(mesh$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = d, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) + 
  labs(x = "Easting (km)", y = "Northing (km)")

# Basic model
m0 <- sdmTMB(
  length_mm ~ temp_sc*species + year_ct*species + yday_ct,
  data = d,
  mesh = mesh,
  family = gengamma(link = "log"),
  spatiotemporal = "off",
  spatial = "on",
  time = "year")

sanity(m0)
m0
Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc * species + year_ct * species + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: d
Family: gengamma(link = 'log')
 
                                        coef.est coef.se
(Intercept)                                 2.18    0.02
temp_sc                                     0.03    0.01
speciesAmmodytidae                          1.69    0.02
speciesAnguilla anguilla                    2.09    0.02
speciesAphia minuta                         1.34    0.01
speciesArgentina silus                      1.58    0.16
speciesArgentina spp                        1.28    0.03
speciesChirolophis  ascanii                 0.60    0.02
speciesClupea harengus                      1.10    0.01
speciesCrystallogobius linearis             0.97    0.01
speciesEnchelyopus cimbrius                 1.04    0.05
speciesLimanda limanda                     -0.39    0.03
speciesMaurolicus muelleri                  1.46    0.04
speciesMicrostomus kitt                     1.05    0.02
speciesMyoxocephalus scorpius              -0.09    0.02
speciesPholis gunnellus                     0.43    0.01
speciesPomatoschistus sp                    1.49    0.01
speciesSardina pilchardus                   1.30    0.08
speciesSprattus sprattus                    1.72    0.02
speciesSyngnathus rostellatus               2.03    0.02
speciesTaurulus bubalis                    -0.22    0.05
year_ct                                     0.00    0.00
yday_ct                                     0.00    0.00
temp_sc:speciesAmmodytidae                  0.02    0.01
temp_sc:speciesAnguilla anguilla           -0.01    0.02
temp_sc:speciesAphia minuta                -0.02    0.01
temp_sc:speciesArgentina silus             -0.04    0.04
temp_sc:speciesArgentina spp               -0.18    0.03
temp_sc:speciesChirolophis  ascanii        -0.02    0.01
temp_sc:speciesClupea harengus             -0.05    0.01
temp_sc:speciesCrystallogobius linearis    -0.03    0.01
temp_sc:speciesEnchelyopus cimbrius         0.19    0.05
temp_sc:speciesLimanda limanda              0.07    0.02
temp_sc:speciesMaurolicus muelleri         -0.09    0.03
temp_sc:speciesMicrostomus kitt            -0.11    0.02
temp_sc:speciesMyoxocephalus scorpius      -0.04    0.02
temp_sc:speciesPholis gunnellus            -0.06    0.01
temp_sc:speciesPomatoschistus sp           -0.01    0.01
temp_sc:speciesSardina pilchardus          -0.10    0.02
temp_sc:speciesSprattus sprattus           -0.03    0.02
temp_sc:speciesSyngnathus rostellatus       0.00    0.01
temp_sc:speciesTaurulus bubalis            -0.15    0.04
speciesAmmodytidae:year_ct                 -0.02    0.00
speciesAnguilla anguilla:year_ct            0.01    0.00
speciesAphia minuta:year_ct                 0.00    0.00
speciesArgentina silus:year_ct             -0.04    0.02
speciesArgentina spp:year_ct               -0.01    0.00
speciesChirolophis  ascanii:year_ct         0.00    0.00
speciesClupea harengus:year_ct              0.00    0.00
speciesCrystallogobius linearis:year_ct     0.00    0.00
speciesEnchelyopus cimbrius:year_ct        -0.01    0.00
speciesLimanda limanda:year_ct              0.04    0.00
speciesMaurolicus muelleri:year_ct          0.02    0.00
speciesMicrostomus kitt:year_ct             0.00    0.00
speciesMyoxocephalus scorpius:year_ct       0.00    0.00
speciesPholis gunnellus:year_ct             0.00    0.00
speciesPomatoschistus sp:year_ct            0.01    0.00
speciesSardina pilchardus:year_ct          -0.01    0.01
speciesSprattus sprattus:year_ct            0.00    0.00
speciesSyngnathus rostellatus:year_ct       0.00    0.00
speciesTaurulus bubalis:year_ct             0.01    0.01

Dispersion parameter: 0.22
Generalized gamma lambda: 0.51
Matérn range: 13.42
Spatial SD: 0.09
ML criterion at convergence: 149690.315

See ?tidy.sdmTMB to extract these values as a data frame.
# IID ST fields
m1 <- sdmTMB(
  length_mm ~ temp_sc*species + year_ct*species + yday_ct,
  data = d,
  mesh = mesh,
  family = gengamma(link = "log"),
  spatiotemporal = "IID",
  spatial = "on",
  time = "year")

sanity(m1)
m1
Spatiotemporal model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc * species + year_ct * species + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: d
Family: gengamma(link = 'log')
 
                                        coef.est coef.se
(Intercept)                                 2.20    0.02
temp_sc                                     0.02    0.01
speciesAmmodytidae                          1.65    0.02
speciesAnguilla anguilla                    2.08    0.02
speciesAphia minuta                         1.31    0.01
speciesArgentina silus                      1.64    0.17
speciesArgentina spp                        1.27    0.03
speciesChirolophis  ascanii                 0.58    0.02
speciesClupea harengus                      1.09    0.01
speciesCrystallogobius linearis             0.96    0.01
speciesEnchelyopus cimbrius                 0.99    0.05
speciesLimanda limanda                     -0.34    0.03
speciesMaurolicus muelleri                  1.43    0.04
speciesMicrostomus kitt                     1.05    0.02
speciesMyoxocephalus scorpius              -0.08    0.02
speciesPholis gunnellus                     0.40    0.01
speciesPomatoschistus sp                    1.47    0.01
speciesSardina pilchardus                   1.26    0.09
speciesSprattus sprattus                    1.70    0.02
speciesSyngnathus rostellatus               2.01    0.02
speciesTaurulus bubalis                    -0.21    0.05
year_ct                                    -0.01    0.00
yday_ct                                     0.00    0.00
temp_sc:speciesAmmodytidae                  0.03    0.01
temp_sc:speciesAnguilla anguilla           -0.01    0.02
temp_sc:speciesAphia minuta                -0.02    0.01
temp_sc:speciesArgentina silus             -0.03    0.04
temp_sc:speciesArgentina spp               -0.17    0.03
temp_sc:speciesChirolophis  ascanii        -0.02    0.01
temp_sc:speciesClupea harengus             -0.05    0.01
temp_sc:speciesCrystallogobius linearis    -0.03    0.01
temp_sc:speciesEnchelyopus cimbrius         0.17    0.05
temp_sc:speciesLimanda limanda              0.08    0.02
temp_sc:speciesMaurolicus muelleri         -0.10    0.03
temp_sc:speciesMicrostomus kitt            -0.12    0.02
temp_sc:speciesMyoxocephalus scorpius      -0.03    0.02
temp_sc:speciesPholis gunnellus            -0.05    0.01
temp_sc:speciesPomatoschistus sp           -0.02    0.01
temp_sc:speciesSardina pilchardus          -0.10    0.02
temp_sc:speciesSprattus sprattus           -0.03    0.02
temp_sc:speciesSyngnathus rostellatus      -0.01    0.01
temp_sc:speciesTaurulus bubalis            -0.12    0.04
speciesAmmodytidae:year_ct                 -0.02    0.00
speciesAnguilla anguilla:year_ct            0.01    0.00
speciesAphia minuta:year_ct                 0.00    0.00
speciesArgentina silus:year_ct             -0.05    0.02
speciesArgentina spp:year_ct               -0.01    0.00
speciesChirolophis  ascanii:year_ct         0.00    0.00
speciesClupea harengus:year_ct              0.00    0.00
speciesCrystallogobius linearis:year_ct     0.00    0.00
speciesEnchelyopus cimbrius:year_ct         0.00    0.00
speciesLimanda limanda:year_ct              0.03    0.00
speciesMaurolicus muelleri:year_ct          0.02    0.00
speciesMicrostomus kitt:year_ct             0.00    0.00
speciesMyoxocephalus scorpius:year_ct       0.00    0.00
speciesPholis gunnellus:year_ct             0.01    0.00
speciesPomatoschistus sp:year_ct            0.01    0.00
speciesSardina pilchardus:year_ct          -0.01    0.01
speciesSprattus sprattus:year_ct            0.00    0.00
speciesSyngnathus rostellatus:year_ct       0.00    0.00
speciesTaurulus bubalis:year_ct             0.00    0.01

Dispersion parameter: 0.21
Generalized gamma lambda: 0.51
Matérn range: 51.48
Spatial SD: 0.04
Spatiotemporal IID SD: 0.09
ML criterion at convergence: 148310.742

See ?tidy.sdmTMB to extract these values as a data frame.

Cross validation

set.seed(99)

clust <- kmeans(d[, c("X", "Y")], 5)$cluster

d$clust <- clust

pal <- coul[c(1, 6, 8, 16, 20)]

plot_map +
  geom_point(data = d, aes(X*1000, Y*1000, color = as.factor(clust)),
             alpha = 0.3, size = 0.6, shape = 21, fill = NA) +
  scale_color_manual(values = pal)
Warning: `guide_colourbar()` needs continuous scales.

ggsave(paste0(home, "/figures/cv_cluster.pdf"), width = 11, height = 11, units = "cm")
Warning: `guide_colourbar()` needs continuous scales.
m0_cv <- sdmTMB_cv(
  length_mm ~ temp_sc*species + year_ct*species + yday_ct,
  data = d,
  mesh = mesh,
  family = gengamma(link = "log"),
  spatiotemporal = "off",
  spatial = "on",
  fold_ids = clust,
  k_folds = length(unique(clust))
)

m1_cv <- sdmTMB_cv(
  length_mm ~ temp_sc*species + year_ct*species + yday_ct,
  data = d,
  mesh = mesh,
  family = gengamma(link = "log"),
  spatiotemporal = "IID",
  spatial = "on",
  time = "year",
  fold_ids = clust,
  k_folds = length(unique(clust))
)

# Compare log-likelihoods -- higher is better!
m0_cv$sum_loglik
[1] -154327.9
m1_cv$sum_loglik
[1] -153179.6

Go with model 1, spatiotemporal IID

Plot residuals

# samps <- sdmTMBextra::predict_mle_mcmc(m1, mcmc_iter = 201, mcmc_warmup = 200)
# d$mcmc_res1 <- as.vector(residuals(m1, type = "mle-mcmc", mcmc_samples = samps))
# #d$mcmc_res1 <- residuals(m1)
# 
# # Looks reasonably, a weird at the tails though
# d %>%
#   ggplot(aes(sample = mcmc_res1)) +
#   stat_qq(size = 0.75, shape = 21, fill = NA) +
#   stat_qq_line() +
#   labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
#   theme(aspect.ratio = 1)
# 
# ggsave(paste0(home, "/figures/qq.pdf"), width = 11, height = 11, units = "cm")
# 
# # Plot in space
# plot_map +
#   geom_point(data = d, aes(X*1000, Y*1000, color = mcmc_res1)) +
#   scale_color_gradient2()

# Residuals -- above doesn't work yet with gengamma!
png(file = paste0(home, "/figures/qq.png"), units = "cm", width = 15, height = 15, res = 300)
simulate(m1, nsim = 300, type = 'mle-mvn') %>%
  dharma_residuals(m1)
dev.off()
quartz_off_screen 
                2

Predict on grid and plot spatial

pred <- predict(m1, newdata = pred_grid %>%
                  mutate(species = "Clupea harengus")) # The ST and S random fields do not change with species

# Spatial random
plot_map +
  geom_raster(data = pred %>% filter(year == max(year)),
              aes(X*1000, Y*1000, fill = omega_s)) +
  scale_fill_gradient2(name = "Spatial random effect") +
  theme_sleek(base_size = 9) +
  theme(legend.position.inside = c(0.25, 0.1),
        legend.direction = "horizontal",
        legend.key.width = unit(0.7, "cm"),
        legend.key.height = unit(0.3, "cm")) +
  guides(fill = guide_colorbar(position = "inside",
                               title.position = "top",
                               title.hjust = 0.5)) + 
  geom_sf()

ggsave(paste0(home, "/figures/omega.pdf"), width = 11, height = 11, units = "cm")

# Spatiotemporal random
plot_map_fc +
  geom_raster(data = pred,
              aes(X*1000, Y*1000, fill = epsilon_st)) +
  scale_fill_gradient2(name = "Spatiotemporal random effects") + 
  facet_wrap(~year, ncol = 7) + 
  theme_facet_map(base_size = 10) + 
  geom_sf()

ggsave(paste0(home, "/figures/epsilon.pdf"), width = 17, height = 15, units = "cm")


# Prediction (herring), selected years
plot_map_fc +
  geom_raster(data = pred %>% filter(year %in% c(1992, 2002, 2012, 2022)),
              aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(name = "Length (mm)", trans = "sqrt") +
  facet_wrap(~year, ncol = 4) +
  theme_sleek(base_size = 9) +
  theme(legend.position.inside = c(0.195, 0.75),
        legend.direction = "vertical",
        legend.key.width = unit(0.2, "cm"),
        legend.key.height = unit(0.3, "cm"),
        legend.title = element_text(size = 6)) +
  guides(fill = guide_colorbar(position = "inside",
                               title.position = "top",
                               title.hjust = 1)) + 
  geom_sf()

ggsave(paste0(home, "/figures/st_pred.pdf"), width = 17, height = 7, units = "cm")

# Plot uncertainty
sim <- predict(m1, nsim = 100, 
               newdata = pred_grid %>%
                 mutate(species = "Clupea harengus"))

sim_last <- sim[pred_grid$year == max(pred_grid$year), ] # just plot last year
pred_last <- pred[pred$year == max(pred_grid$year), ]
pred_last$lwr <- apply(exp(sim_last), 1, quantile, probs = 0.025)
pred_last$upr <- apply(exp(sim_last), 1, quantile, probs = 0.975)
pred_last$sd <- round(apply(exp(sim_last), 1, function(x) sd(x)), 2)
pred_last$cv <- round(apply(exp(sim_last), 1, function(x) sd(x) / mean(x)), 2)

plot_map +
  geom_raster(data = pred_last,
              aes(X*1000, Y*1000, fill = cv)) +
  scale_fill_viridis(option = "mako") + 
  theme_sleek(base_size = 9) +
  theme(legend.position.inside = c(0.25, 0.1),
        legend.direction = "horizontal",
        legend.key.width = unit(0.7, "cm"),
        legend.key.height = unit(0.3, "cm")) +
  guides(fill = guide_colorbar(position = "inside",
                               title.position = "top",
                               title.hjust = 0.5)) + 
  geom_sf()

ggsave(paste0(home, "/figures/spatial_cv.pdf"), width = 11, height = 11, units = "cm")

Coefficients

tidy(m1) %>% as.data.frame()
                                      term      estimate    std.error
1                              (Intercept)  2.2017933324 0.0170624303
2                                  temp_sc  0.0246658086 0.0128297468
3                       speciesAmmodytidae  1.6468346729 0.0160873369
4                 speciesAnguilla anguilla  2.0773744737 0.0194810177
5                      speciesAphia minuta  1.3087453316 0.0136375644
6                   speciesArgentina silus  1.6408805276 0.1678850380
7                     speciesArgentina spp  1.2661382522 0.0303506268
8              speciesChirolophis  ascanii  0.5844242365 0.0157352162
9                   speciesClupea harengus  1.0942305689 0.0135770762
10         speciesCrystallogobius linearis  0.9636805608 0.0136184992
11             speciesEnchelyopus cimbrius  0.9921311106 0.0469485503
12                  speciesLimanda limanda -0.3413526086 0.0307652385
13              speciesMaurolicus muelleri  1.4316224155 0.0400275537
14                 speciesMicrostomus kitt  1.0468385677 0.0165714706
15           speciesMyoxocephalus scorpius -0.0795154997 0.0191302853
16                 speciesPholis gunnellus  0.4031977450 0.0136823916
17                speciesPomatoschistus sp  1.4651207322 0.0142411280
18               speciesSardina pilchardus  1.2560352585 0.0909112746
19                speciesSprattus sprattus  1.6979418847 0.0176364136
20           speciesSyngnathus rostellatus  2.0109400867 0.0157665100
21                 speciesTaurulus bubalis -0.2103747533 0.0521629182
22                                 year_ct -0.0051828622 0.0013064920
23                                 yday_ct  0.0018833563 0.0005274263
24              temp_sc:speciesAmmodytidae  0.0312250113 0.0138272361
25        temp_sc:speciesAnguilla anguilla -0.0072353640 0.0202931399
26             temp_sc:speciesAphia minuta -0.0204725003 0.0121542931
27          temp_sc:speciesArgentina silus -0.0258507957 0.0374055255
28            temp_sc:speciesArgentina spp -0.1744859525 0.0335999122
29     temp_sc:speciesChirolophis  ascanii -0.0226885606 0.0136453105
30          temp_sc:speciesClupea harengus -0.0493219672 0.0123029235
31 temp_sc:speciesCrystallogobius linearis -0.0334516293 0.0121512676
32     temp_sc:speciesEnchelyopus cimbrius  0.1705620970 0.0459767541
33          temp_sc:speciesLimanda limanda  0.0800843872 0.0184944242
34      temp_sc:speciesMaurolicus muelleri -0.0958958240 0.0316982103
35         temp_sc:speciesMicrostomus kitt -0.1187109721 0.0168959849
36   temp_sc:speciesMyoxocephalus scorpius -0.0270361945 0.0152627963
37         temp_sc:speciesPholis gunnellus -0.0534576875 0.0120790401
38        temp_sc:speciesPomatoschistus sp -0.0167895195 0.0125776493
39       temp_sc:speciesSardina pilchardus -0.1020761875 0.0190106294
40        temp_sc:speciesSprattus sprattus -0.0323183759 0.0155200946
41   temp_sc:speciesSyngnathus rostellatus -0.0137953804 0.0137118139
42         temp_sc:speciesTaurulus bubalis -0.1201217321 0.0429687448
43              speciesAmmodytidae:year_ct -0.0169662389 0.0014053404
44        speciesAnguilla anguilla:year_ct  0.0053383930 0.0019758344
45             speciesAphia minuta:year_ct  0.0015989702 0.0012125807
46          speciesArgentina silus:year_ct -0.0497146976 0.0176120858
47            speciesArgentina spp:year_ct -0.0079958455 0.0022274962
48     speciesChirolophis  ascanii:year_ct -0.0012728807 0.0013582374
49          speciesClupea harengus:year_ct  0.0015120757 0.0011879410
50 speciesCrystallogobius linearis:year_ct -0.0002908017 0.0011894349
51     speciesEnchelyopus cimbrius:year_ct -0.0013727965 0.0034116001
52          speciesLimanda limanda:year_ct  0.0320642189 0.0027933830
53      speciesMaurolicus muelleri:year_ct  0.0246322658 0.0035843246
54         speciesMicrostomus kitt:year_ct -0.0009833155 0.0015857900
55   speciesMyoxocephalus scorpius:year_ct -0.0002173287 0.0018444513
56         speciesPholis gunnellus:year_ct  0.0054458029 0.0012005775
57        speciesPomatoschistus sp:year_ct  0.0077844893 0.0012831516
58       speciesSardina pilchardus:year_ct -0.0090907467 0.0059110581
59        speciesSprattus sprattus:year_ct -0.0024785831 0.0017975039
60   speciesSyngnathus rostellatus:year_ct  0.0006984562 0.0015935258
61         speciesTaurulus bubalis:year_ct  0.0042812992 0.0087080424
sort(unique(d$species))
 [1] "Agonus cataphractus"      "Ammodytidae"             
 [3] "Anguilla anguilla"        "Aphia minuta"            
 [5] "Argentina silus"          "Argentina spp"           
 [7] "Chirolophis  ascanii"     "Clupea harengus"         
 [9] "Crystallogobius linearis" "Enchelyopus cimbrius"    
[11] "Limanda limanda"          "Maurolicus muelleri"     
[13] "Microstomus kitt"         "Myoxocephalus scorpius"  
[15] "Pholis gunnellus"         "Pomatoschistus sp"       
[17] "Sardina pilchardus"       "Sprattus sprattus"       
[19] "Syngnathus rostellatus"   "Taurulus bubalis"        
mean_coefs <- tidy(m1, conf.int = TRUE) %>%
  filter(term %in% c("temp_sc", "yday_ct")) %>% 
  mutate(term = ifelse(term == "temp_sc",
                       "temp_sc:speciesAgonus cataphractus",
                       "speciesAgonus cataphractus:year_ct"))

# Wrangle data
coefs_temp <- tidy(m1, conf.int = TRUE) %>%
  filter(grepl("temp_sc:", term)) %>% 
  bind_rows(mean_coefs %>% filter(term == "temp_sc:speciesAgonus cataphractus")) %>% 
  mutate(estimate2 = mean_coefs %>%
           filter(term == "temp_sc:speciesAgonus cataphractus") %>%
           pull(estimate) + estimate,
         
         conf.low2 = mean_coefs %>%
           filter(term == "temp_sc:speciesAgonus cataphractus") %>%
           pull(conf.low) + conf.low,
         
         conf.high2 = mean_coefs %>%
           filter(term == "temp_sc:speciesAgonus cataphractus") %>%
           pull(conf.high) + conf.high) %>% 
  mutate(species = str_remove(term, "temp_sc:"),
         species = str_remove(species, "species"))


coefs_year <- tidy(m1, conf.int = TRUE) %>%
  filter(grepl(":year_ct", term)) %>% 
  bind_rows(mean_coefs %>% filter(term == "speciesAgonus cataphractus:year_ct")) %>% 
  mutate(estimate2 = mean_coefs %>%
           filter(term == "speciesAgonus cataphractus:year_ct") %>%
           pull(estimate) + estimate,
         
         conf.low2 = mean_coefs %>%
           filter(term == "speciesAgonus cataphractus:year_ct") %>%
           pull(conf.low) + conf.low,
         
         conf.high2 = mean_coefs %>%
           filter(term == "speciesAgonus cataphractus:year_ct") %>%
           pull(conf.high) + conf.high) %>% 
  mutate(species = str_remove(term, ":year_ct"),
         species = str_remove(species, "species"))

p1 <- coefs_temp %>% 
  mutate(sign = ifelse(estimate2 > 0, "pos", "neg"),
         sig = ifelse(estimate2 > 0 & conf.low2 > 0, "sig", "not sig"),
         sig = ifelse(estimate2 < 0 & conf.high2 < 0, "sig", sig)) %>% 
  ggplot(aes(estimate2, reorder(species, coefs_year$estimate2), fill = sig, color = sign, shape = sig)) + 
  geom_point(fill = NA) + 
  scale_shape_manual(values = c(21, 19)) +
  geom_errorbar(aes(xmin = conf.low2, xmax = conf.high2), 
                width = 0, alpha = 0.3) + 
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2, linewidth = 0.25) +
  theme(axis.text.y = element_text(face = "italic")) + 
  labs(y = "Species", x = "Temperature slope") + 
  scale_color_brewer(palette = "Dark2", direction = -1) +
  scale_fill_brewer(palette = "Dark2", direction = -1) +
  guides(color = "none", shape = "none") +
  geom_stripped_rows(aes(y = species), inherit.aes = FALSE) +
  theme(legend.position = "bottom")

p2 <- coefs_year %>% 
  mutate(sign = ifelse(estimate2 > 0, "pos", "neg"),
         sig = ifelse(estimate2 > 0 & conf.low2 > 0, "sig", "not sig"),
         sig = ifelse(estimate2 < 0 & conf.high2 < 0, "sig", sig)) %>% 
  ggplot(aes(estimate2, reorder(species, coefs_year$estimate2), fill = sig, color = sign, shape = sig)) + 
  geom_point(fill = NA) + 
  scale_shape_manual(values = c(21, 19)) +
  geom_errorbar(aes(xmin = conf.low2, xmax = conf.high2), 
                width = 0, alpha = 0.3) + 
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2, linewidth = 0.25) +
  theme(axis.text.y = element_text(face = "italic")) + 
  labs(y = "Species", x = "Year slope") + 
  scale_color_brewer(palette = "Dark2", direction = -1) +
  scale_fill_brewer(palette = "Dark2", direction = -1) +
  guides(color = "none", shape = "none") +
  geom_stripped_rows(aes(y = species), inherit.aes = FALSE) +
  theme(legend.position = "bottom")

p2 + p1 + plot_layout(axes = "collect")

ggsave(paste0(home, "/figures/coefs.pdf"), width = 17, height = 10, units = "cm")